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ABSTRACT 


A computer simulation study of channeling in a diamond 
lattice. The simulation was done for a xenon ion striking 
Ene (100) 4110) or (111) surtace ona silicon tarcaenw. 
Potential functions for the Si-Si lattice bond and the Xe-Si 
interaction are postulated. The electronic stopping cross 
section for the (110) channel of silicon is estimated. 

This work is a continuation in the development of a 
computer model formulated at the USNPGS which takes into 
consideration the displacement of the atoms in the target 
lattice as well as inelastic energy losses by the primary 
ion. The lattice was not thermalized and only the repulsive 
portion of the lattice-lattice potential was utilized. 


Computed ranges are in good agreement with experimental date. 
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Range distribution curves for 40 Kev Xe in 
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Radial Electron Density functions for a number of 
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IPF's for various identical ion interactions of 
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for comparison. 

Radial Electron Density functions for the neutral 
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Range of a 5 KevV Xe 10on at the various 
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I. INTRODUCTION 


Experimental and theoretical studies over the past decade 
have established that the slowing down of ions in crystals is 
strongly dependent on the crystallographic orientation of the 
lattice with respect to the direction of the incident beam. 
(1,2,3,4,5] This effect is known as channeling. A particle 
moving along one of the major axis of a crystal may be "steered! 
(i.e. its momentum redirected toward the center of the channel) 
by successive gentle collisions with the "strings" of atoms that 
comprise the "walls'' of the channel. Davies et al [1] reported 
that this channeling effect increases the mean range of wena 
by as much as a factor of 10 over the values observed ina 
non-crystalline solid of similar atomic number. 

Earlier computer studies performed at the NPGS investigated 
channeling in the body centered [6] and face centered cubic 
ihattace. [7] This thesis is an extension of this previous work 
imtoo the diamond laittice. 

When a substrate is bombarded by a beam of incident ener- 
getic particles it will lose some of its own ions by sputtering 
and retain some of the incident ions. The incident ions retained 
are said to be implanted, and the technique of using an energetic 
10n beam to introduce ions into a substrate is called ion implan- 
potion. | Sie. This doping technique in silicon should be very 
effective due to the extremely open nature of its lattice and 


its large channel Sizes. 


2) 


The difference in the depth of penetration of energetic 
particles along the prinedpal@channels in silicon mayebe attra 
buted to the various channel sizes. Figures 1-3 show the loca- 
tions of the various channels in the (100), (110) and (111) 
surfaces of silicon while Fig. 4 is a "ball and stick" model 
showing the relative shape and sizes of the three principal 
channels in silicon. According to Davies [1] and Manchester [8] 
the preferential order of channeling would occur along the (110), 
(111) and (100) channels in decreasing order with the (111) 
channel only marginally better than the (100) channel. This fact 
is graphically indicated in Fig. 5. 

Much work over the past several years has been devoted to 
studying this channeling effect in silicon. Manchester | 8] has 
Studied the feasibility of forming metallurgical junctions in 
Silicon by making use of the channeling phenomenon. Eisen|9} 
has studied the channeling of medium mass ions through silicon 
with special attention directed toward the electronic stopping 
cross section of well channeled ions. Davies and his co-workers 
have written a series of papers [2,3,4,10] discussing the pene- 
tration of KeV and MeV projectiles in silicon and the disordering 
of the lattice caused by this ion bombardment. Glotin Cae 
Nelson and Mazey hee and Dearnaley et al [13] have investigated 
the implantation of phosphorous in silicon. Eriksson et al [14] 
have studied the implantation and annealing behavior of Group III 
and V ions in silicon. Gibbons [5] sums up the current theoreti- 
cal and experimental work pertinent to the problem of predicting 


impurity distribution profiles in implanted material. 
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Ii. STUDY OBJECTIVES 


This report is an extension of the NPGS channeling model from 
the bec and fee orientations into the diamond lattice. This 
thesis attempts to establish, with xenon as the bullet and sili- 
con as the target, the following characteristics for the three 
major channeling orientations: 

a) The potential between the bullet and the lattice, 

b) A value for the electronic stopping power, and 


c) The ion implantation profile at 5 KeV. 
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III. SIMULATION MODEL 


The basic model has been developed by the NPGS group over 
the past six years. A brief description of this model and the 
alterations made for this study follow but for detailed pre- 
sentations the reader is referred to Refs.7,15 and 16. 

The simulation model consists of a single primary anda 
silicon (diamond) lattice target containing between 60 and 150 
atoms. All runs were made with the primary (bullet) striking 
the target normally, although the program is sufficiently 
flexible that the bullet could be fired into the target at any 
angle. 

The silicon target lattice as used in this simulation is 
composed of a central core and a surrounding shell. The program 
is structured so that the bullet is never allowed to leave this 
central core. In addition, in the (110) and (111) orientations 
the shell atoms were not permitted to move. This allows the 
program to run slightly faster and also tends to hold the core 
atoms in place for a slightly longer time. The shell atoms were 
not allowed to displace in the (100) channeling direction. Even 
with the shell atoms displacing, given a bullet energy and an 
equivalent channeling location in the (111) and (100) orien- 
tations, the (100) channeling program ran faster. For a detailed 
discussion of the lattice generator to inelude the dividingyor 
the target crystal into active core and passive shell see 


Appendix A. 
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In previous computer simulations the interatomic potential 
function and force function were of the exponential (Born-Mayer ) 
type 


F = exp(A + Bx) 


where A and B are empirically determined constants and x 1s the 
atomic separation. The functions are not always of this type 
in this simulation, and force and potential tables were con- 
structed. For details of this scheme see Appendix B. 

Two major assumptions are made in this simulation model, 

(1) The bullet moves in a perfect lattice undisturbed by the 
thermal displacements of the atoms. 

(2) Only repulsive forces between the lattice atoms have 
been included. 

While both of these assumptions are obviously erroneous, 
previous studies have shown that this simulation model is still 


a good first order approximation for channeling. 


is 


IV.” PROCEDURE 


This simulation study was conducted in an IBM 360/67 computer 
using FORTRAN IV language. 

The computer sets up the diamond lattice in the desired 
orientation (see Figs. 1-4) and starts a xenon ion into it. A 
brief study of the symmetry of the target lattice shows that the 
indicated impact areas cover all the possible points for the 
three orientations of silicon. Within these impact areas, 
several points, depending on the orientation of the lattice, were 
selected as representative of the channel. In the determination 
of how many points were to be chosen in each impact area, con- 
Sideration was given to obvious lattice symmetry and to computer 
running time. 

Various combinations of Bullet-Target Potentials and CELS, 

a constant used to determine inelastic energy losses, were tried 
at different bullet energies and the maximum range of the bullet 
for each orientation was matched against experimental range data 
as reported by Davies et al [1]. 

Channel profiles for 5 KeV xenon atoms were determined for 
‘the (100), (110), and (111) channels of silicon. Runs were made 
to determine dE/dx for each major axis in silicon and maximum 
ranges for xenon at higher energies were determined knowing 
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V. POTENTIAL FUNCTIONS AND INELASTIC LOSS CONSTANT 


The purpose of this section is to discuss first, the 
determination of the potentials used to approximate the Si-Si 
lattice bond and the xenon-silicon bullet-lattice interaction, 
and second, the determination of an inelastic loss constant. 
Presently, there exists in the literature calculated potentials 
for the Si-Si system [17] but none for the Xe-Si system. 

Eisen [9] has published experimental electronic stopping cross 
sections for medium mass ions in silicon but no data exists for 


xenon in silicon. 


A. S$1i-Si BOND 

The starting point for selecting a potential for the S$i-Si 
bond are the Interaction Potential Functions (IPF) as calculated 
by Harrison [17]. These potentials are only the repulsive portion 
of the potential function for pairs of particles and neglect any 
lattice effects. Harrison [17] also indicates that they are most 
likely too "hard". Figure 6 shows the radial electron density 
for Si. Since this radial electron density was determined using 
a spherical approximation, it is not correct for r 2 1.5a,, 
(a, = Bohr Radius = 0.5294) Eis Figure 7 shows several dif- 
ferent potential functions for the Si-Si bond with the Gibson II 
sacs comparison. The Si-Si and Si-Si. IPF's show approximately 
the correct behavior and zero close to the nearest neighbor dis- 


—: : As . 
tance for silicon. The Sr si : also shows approximately the 
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correct behavior but fails to zero at the nearest neighbor 
. . ; carl +1 
distance. Curve 3 of Fig. 7 1s the Si -S1 IPF forced to 
zero at the nearest neighbor distance. 
2 . : , Pa f oe 
Any one of these three potentials (S1-S1, Si-Si or Si-Si 
(zeroed)) will serve as a good approximation for the Si-Si bond. 
It is obvious that the Gibson II will not serve as a good repre- 
sentation of this bond. 
The covalent bonding between Silicon atoms has the effect 
of concentrating about 0.5 of an electron of the outer shell 
along every chemical bond iegeger Therefore, it seems reasonable 
to conclude that each silicon atom along the Si-Si bond sees the 
other in a slightly ionized state but not completely ionized 
to the +l state due to the shielding effects of this shared 
: a : : +1 one 
electron. This would seem to eliminate the S1-S1 and S1 §$-S1 
potentials for the Si-Si bond and lead to the selection of the 
+ ; : 
Si-Si potential as representing the bond. It was noted earlier, 
however, that the IPF's are probably too hard. Therefore, the 
Weeli +1 : : 
Si. 251 (zeroed) potential was chosen as representative of the 
: or ; 
bond potential as it is "softer" than the Si-Si potential and 
is forced to zero at the nearest neighbor distance. 
= + + 
Admittedly the reasons for choosing the S1 -S1 potential 
were rather arbitrary but the reader 1s asked to keep three 
points in mind. First, all three potentials are approximately 
equal and trial runs indicated no significant effect on the range 
when these potentials were interchanged. Second, the neutral 


Silicon atoms,are toomlarge topfiateinto the Sisbattaee ahd, 
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the two overwhelming factors that control channeling ranges are 
the bullet-lattice potential and the inelastic loss approximation. 
Any reasonably potential selected for the Si-Si bond should have 


no discernable effect on these ranges. 


B. INELASTIC LOSS CONSTANT 

One of the principal mechanisms of energy loss for channeled 
pons im a sold gs the interaction of the Yvon with electrons 
(both bound and free) of the solid. This computer study accounts 
for this loss through the use of a "frictional force multiplier" 
which acts on the ions velocity in a manner to produce a "drag"! 
merce On the incident i10n. Thas "fractaenefonce.multaphzer™ is 
directly proportional to the electronic stopping cross section. 

Eisen [9] reports electronic stopping cross sections (S.) 
f6r medium mass ions (5 = Z < 19) in the (100), (110), ana (114) 
channels of silicon. In unpublished data [19] he extends Ss. for 
the (110) channel by allowing the atomic number of the incident 
ion to increase up to Z = 36 (Krypton). Eisen's data reveals the 
felblowang: first, Ss. was found to exhibit a strong oscillatory 
dependence on Z, the atomic number of the incident ion, second, Ss. 
is lowest in the (110) channel and third, S, is approximately 
equal in the ¢111) and (100) channels. 

No data presently exists in the literature for the electronic 
stopping cross section of xenon in silicon. However, based on 
previous computer studies of channeling in copper and tungsten 
undertaken at NPGS the frictional force multiplier was estimated 


fomeewOs10772 hesec/m in te {100) channel. Using Eisen's 


ie 


data [19] and noting the ratio between S. in the (100) and (110) 
channels, the frictional force multiplier in the €110) channel 


Jia 


is approximately -8 x 10° n-sec/m. These values were used in 
this computer study. 

Ericksson et al [20] give estimates of electronic stopping 
cross eeoten for various ions along the (100) and (110) channels 
of tungsten. Their data also indicates that S. shows a strong 
oscillatory dependence on Z, the atomic number of the incident 
lon and.-this oscillatory dependence is very similar to the 
dependence Eisen [9,19] found for ions channeled in silicon. If 
the assumption is made that the ratio between stopping cross 
sections for xenon and krypton (which is chemically similar to 
xenon) channeled in tungsten is approximately the same as the 
ratio between these ions channeled in silicon, one obtains a 


value of 


oy 
= aevy em 


8 
rn (where v = 1.0 x 10° cm/Sec) 


8.5 x 10 


as an estimate of Ss. inet hes iMl@eie haname is lnmesrniie onpmchh 1s 
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converts to a value of 7 x 1:Ow. n-sec/m as a frictional foxide 
multiplier for the (110) channel and is in good agreement with 


the value used in this computer study. 


S aye -19 
CELS = —2* ea er esee oaeN where 7] = number density. 


Sm eran 


v(cm/sec) x 10. 


The author does not make any claim that the values of the 
“Frictional force!smused anpthisesctudyeare “corr come Tt 1s fear. 


however, that based on a study of Eisen's data ROR19) and the 
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extrapolation described above, the values chosen for S. will 
certainly be within 20% of the experimentally determined values 


and more than likely will be much closer than that. 


C. Xe-Si POTENTIAL 

As in the case of the Si-Si bond, no potential function 
exists for the xenon-Silicon interaction. Therefore, the 
Interaction Potential Functions will again serve as a starting 
point for selecting the Xe-Si potential. Figure 8 shows the 
radial electron density for xenon in the nonionized and +6 


+l 


Iie +o .. 
1Onization state, while Fig. 9 shows the IPF for Xe -Si1 and 


Far from the S1i-Si bond, the electron density for silicon 


fh 


Should approach the Si’ distribution (Fig. 6). Special note 
should be taken of the €111) channel of silicon where the atoms 
bounding the channel have one of their bonds parallel to the 
centerline of the channel. Determing the ionization state of 
xenon is also difficult. Harrison [17] has found for Xe-W 
system a good potential fit can be found if xenon is assumed to 
be in the +6 or +8 ionization state, even though this high ionic 
state 1S somewhat surprising. The ot i potential was found 
to be unrealistically strong. As a result, of the above it was 


+h 


a , : 
concluded that the Xe -S1 IPF was a good starting point for 


the bullet-lattice potential. 


+6 + 
The Xe ~-S1 4 


IPF (Fig. 9) diverges to unrealistically 
Strong values at large separations. This divergence is charac- 


teristic of the theoretical model used to calculate the IPF 


in 


and is not characteristic of the physical system | 1713 eA Boerne 
Mayer approximation for the Xe-Si interaction would fail at small 
Separations because it would approach a finite value at zero 
separation, while the actual Xe-Si function should approach a 
pure coulomb potential at small separations. Our final choice 

of the potential to represent the Xe-Si system incorporated both 
of the above approximations. For separations less than 1.2a,, 
the IPF was used, while a Born-Mayer approximation was constructed 
for separations greater’ than 1.2a,.. 

Curve BM1l (Fig. 9) when used in conjunction with the pre- 
viously discussed values of the "frictional force" parameter 
gives excellent results for xenon ions injected in the (110) and 
(100) channels, however, it is a poor approximation for xenon 
in the {11l) channel. After several trial runs it became obvious 
t@rat all thee channeling diWections could nowt Prt on “one™@™5een-— 
Mayer approximation. The maximum range in the (111) channel 
could be forced to approach the experimental range if the 
"frictional force" constant is decreased by two orders of magni- 
tude. this 1S UnrealisStic since it conflicts drastically wren 
Eisen's [9] experimental data. 

A solution to the problem was to make a separate Born-Mayer 
approximation for the (111) channel; this is curve BM2 in Fig. 9. 
Intuitively, one feels that this answer is not the correct one. 
All previous studies have found that one potential satisfies all 
channeling directions. However, two points should be kept in 
mind. First, all previous computer studies dealt with channeling 


in a metallic target where the principal method of bonding between 
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lattice atoms is metallic. This is the first computer study where 
the bonding between the lattice atoms is co-valent and any con- 
clusions drawn for a metallic target need not necessarily be true 
in silicon. Secondly, the orientation of the bond between silicon 
atoms with relation to the channeling directions varies among all 
three principal channels of silicon. As previously mentioned, 
the extreme is reached in the (111) channel where the Si-Si bond 
1s parallel to the center line of the channel. This more than 
likely has an effect on the silicon charge distribution as seen 
by the xenon ion and could possibly cause the potential for the 
(111) channel to differ from the other channels. 

Another possibility is that a Born-Mayer approximation is not 


+ : 
4 potential might be 


eee =O. ae 
valid in silicon. A picture of the Xe -Si 
Obtained by combining the portion of the three curves: IPF, BM1l 
and BM2 that we know give us accurate results. Such a composite 


curve is shown in Fig. 10. It is obvious that this composite curve 


is unrealistic. 


Fal 


Lee RESULTS 


A. {1005 CHANNEL 

Proure ll is@a@pilot of dE/dx( in’ /Lu) versus 10n energy for 
a perfectly channeled xenon ion in the (100) direction. (Note 
ILU = 2.714 = Guaanetag/en- Between 30-100 KeV, dE/dx is 
approximately proportional to the energy indicating that in this 
region the elastic and inelastic collision processes are both 
exerting a significant influence. This agrees with experimental 
data reported by Davies [1]. 

Figure 12 is a plot of dE/dx versus B2, Above 800 Kev 
dE/dx™is proper t#onal “ro 5? indicating that as expected the in- 
elastic process is the primary loss mechanism. Both Figs. 11 
and 12 show that below 20 KeV the elastic collision process is 
primarily responsible for slowing down the i0n. 

Figures 13 and 14 are range profiles and range profile 
sections for the (100) channel. Figure 13 shows only one-half 
of the channel, the range profile for the other half is very 
camilate but not Mdenticaletce= this, 

Table 1 (below) is a comparison between experimental ranges 
obtained by Davies [1] and computer calculated ranges for the 
(100) channel. The column labeled dE/dx is the range value 


obtained by using the dE/dx curve in Fig. 11 and the relationship 
E 
-1 
R al ( dE/dx ) dE. 
O 


ae 


The computer calculated range at 40 KeV, although higher than 
the eperimental value reported by Davies, falls within his +3% 


error estimate. 


Range (100) 


Energy Exper imental dE/dx Computer 
5 -- 60 69.7 

20 491.7 485 490.2 
LO G27 1230 1193.8 


Table I. Energy in KeV, range in LU 


*Experimental range for 20 KeV extrapolated 
from Davies ( 4) experimental data. 


In a bce or fcc crystal an 10n injected into the center of 
the channel remains there. This is not the case for an ion 
mepected into the center of the (100) chammel of siiiicon. THe 
10n "wobbles" slightly around the channel center oscillating 
between (110) planes. (Fig. 15) Note the similarity of the 
orbit in the second and fourth quadrant. 

Figure 6 is a Range Profile Section for a 5 KeV silicon ion 
injected into the (100) channel of silicon. Although no experi- 
mental data exists to compare maximum ranges it is interesting to 
note that the range profile section 1s very similar to the range 


profile of xenon in silicon. 


B. (110) CHANNEL 
Figure 7 1s a plot of dE/dx versus energy at various values 


of CELS (analogous to electronic stopping power) for a xenon 10n 
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injected into the center of the (110) channel. It is apparent 
that the energy loss 1S very sensitive to CELS especially at 
higher energies indicating that electronic stopping is the 
principle mechanism of energy loss for well channeled ions in 
the (110) direction. Figure 1571s amplot of dE/dx versus 
alors ee for well channeled ions with CELS = - 8 x 10714 
n-sec/m, the value selected for use in this study. dE/dx is 
proportional to Be at incident ion energies greater than 2 KeV 
again indicating that electronic stopping is the principle loss 
mechanism. 

Figure 19 and Fig. 20 are a range profile and a range profile 
section for the (110) channel. Both show only half the channel 
although the other half of the channel should be identical. 

While it is apparent that the (110) channel is a preferred 
channeling direction the maximum range falls off rapidly if the 
i0n is injected slightly off center. (For example, if the ion 
1s injected approximately O0.2LU from the center of the channel 
the maximum range decreased by one-half). Of particular interest 
is point 1 in the Range Profile Section (Fig. 20). An ion in- 
jected at this point in the channel has a range of 100LU which 
is much greater than expected. The xenon ion injected here is 
strongly influenced by the silicon atom located at position C 
which forces it into the center of the channel and adds signi- 
Picatiuly tOrits range. 


Table II is again a comparison between ranges obtained 


experimentally by Davies [1] and computer calculated ranges. 
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Range 110 


Ener gy Exper imental dE/dx Computer 
5 == 1090 1094 
20 2371 2400 ao 
LO 3812 3800 _ 


Table II. Ranges in LU, Energies in KeV. 


The ranges calculated using dE/dx are in good agreement with 
experimental ranges obtained by Davies [1]. No attempt was made 
to channel a 20 KeV Xenon ion as it would be too expensive in 
computer running time (approximately 90 minutes). 

Figure 21 shows, the ranges obtained for selected points 
in the (110) channel when the energy of the incident xenon ion 
was increased to 20 KeV. Two points should be noted, First, 
the channel "expands" at higher energy as expected. Second, a 
comparison of the range of the ion incident on point 1 with the 
range obtained when the incident ion energy was 5 KeV (Fig. 20) 
Shows that the large range for the 5 KeV ion is a low energy 


phenomenon. 


C. €111) CHANNEL 

Figure 22 is an enlargement of the (111) channel. Initially, 
point 1 was thought to be the center of the channel but subse- 
quent probing of the channel proved this assumption to be incor- 
‘rect. The points numbered 2 and 3 are actually the channel 
"centers". If atoms A, B, Dand B, C, D are thought to form 


triangles, as indicated, then points 2 and 3 may be considered 
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the center of their respective triangle equidistant from the atoms 
that form the triangle's vertices. Additionally, although the 
basic shape of the channel remains unchanged, the atoms that 

bound the channel do not necessarily lie in the planes indicated 
in Fig. 22. Two other possible configurations are: (a) atoms A 
and C in the front plane, B in the third plane and D in the fifem 
plane or (b) atoms A and C in the fifth plane back, D in the third 
plane back and B in the front plane. These three channel con- 
figurations are shown in Fig. 3. They affect the ion only as it 
frrviamiyewernters™the lattice, by determining in which direction 

it will be initially steered and have almost no effect (5 0.5%) 

on the calculated maximum range. 


Figure 23 and 24 are plots of dE/dx vs ion energy, E and 


wer 


BE Up to 100 KeV both elastic and inelastic collisions are 
both important loss mechanisms. Below 20 KeV, elastic losses 
appear to predominate. 

Figure 25 shows dx/dE curves for xenon in the (111) channel. 
The area under these curves should be an accurate representation 
of the total range travelled by an incident xenon ion. The solid 
curve (I) was obtained by injecting xenon ions into the center of 
the (111) channel at various energies, allowing them to penetrate 
the lattice for a short distance (S 5LU) and determining the 
energy loss of the xenon ion per lattice unit. The dotted curve 
(II) iS actually a 20 KeV run for a xenon ion initially injected 
in the center of the (111) channel and followed until it finally 


stops in the crystal. The difference between these two curves 


can be accounted for by noting that an ion injected into the 
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center of the (111) channel deviates considerably from the 
channel center (up to .O9LU) as it passes thru the crystal and 
that this deviation increases as the ion energy decreases. 
Therefore, the implicit assumption made in determining curve I, 
that the ion remains in the center of the (111) channel as it 
passes through the crystal is incorrect and will introduce large 
errors (up to 20%) in ranges at the low end of the energy spec- 
trum where the elastic loss mechanism predominates. Table III 
is a summary of maximum ranges obtained at various energies for 


the (111) channel. 


Range (111) 


Energy Exper imental Curve I Curve II Computer 
5 _ 37 60 D3 
20 511 417 500 511 
4O LZ@9 1050 1170 -- 


Table III. Energy in KeV, range in LU. 


Figure 26 is a representation of the (111) channel and shows 
ranges obtained at different initial impact points in the channel. 
Two things can be noted about the range distribution. First, 
there is a definite peak at the center of the channel and 
secondly, the total range is extremely dependent on exact im- 
pact point in the channel. This last fact makes it impossible 


to draw any meaningful range profiles or ranges profile section. 
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VIL. CONCLUSIONS 


1. Maximum channeled ranges were fit to experimental data 
at 20 KeV. 

2. For the (100) and (111) channels both elastic and 
inelastic losses are important between 20 and 80 Ke\V. 

3. Inelastic losses become predominant at energies greater 
than 400 KeV for the (100) channel. 

4. Inelastic loss seems to be the predominant loss mechanism 
for well channeled ions in the (110) channel even at low energies 
(~2 Kev). 

hb 4 86Thé®ram@e in the GM’ channels @Mtremely sensitive to 


impact parameter and appears to possess no symmetry. 


-14 ne 


See a) Wave (110) 


Ga S was determined to be 8.5 x 10 
Ehoanmeimor silicon. in thesGloe). and <1Lll)wchaunels 


2 
-14 eV-cm 
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7. All three channeling directions for xenon in silicon can 
not be determined using the same Born-Mayer Potential Function. 
A minimum of two Born-Mayer Potentials are required. 

8. At low energies the (100) channel is the preferential 
direction for channeling when compared to the {111) channel. 
This is the reverse of the situation for incident ions energies 


greater than 20 KeV. 
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9. Although the potentials and the "frictional force" con- 
stant were discussed separately it can not be overemphasized 
that all three of these variables must be treated as an inte- 
grated whole. A change in either the bullet-target potential 
or the frictional force constant has a profound effect on the 
maximum channeled range especially in the (110) orientations. 
The final selection of the Born-Mayer approximations to the 
Xe-Si interactions, the frictional force parameter and S$i-S1 


bond potential seems to represent the best of all the parameters. 
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APPENDIX A (LATTICE GENERATOR) 


Previous computer studies undertaken by the NPGS Group over 
the past several years have dealt exclusively with the body 
centered and face centered cubic lattices. In building the 
diamond lattice, which is basically two interpenetrating face 
centered cubic lattices, the primary aim was to construct the 
lattice in such a way that only mimimum modifications on the 
basic computer model for channeling would be necessary. The 
diamond lattice generator was constructed so that no major change 
in the basic lattice regenerator was necessary. Therefore, to 
understand the lattice generator, it is necessary to understand 
the fundamental mechanism behind the regenerator subroutine. 

Subroutine REGEN works as follows: 

(1) It receives the parameter MSS from the main program 
which indicates the direction of regeneration (+X, +Y, +Z). 

(2) It checks every atom in the crystal against a constant 
which depends on the direction of regeneration and either 

a. discards the atom since the bullet has passed and is no 
Longer under its influence, or 

b. shifts and renumbers the atom while preserving its 
velocity components. 

(3) It builds and undisturbed section crystal in the 
direction of travel of the bullet. 

(4) It shifts the bullet to its proper location in the 


rebuilt crystal preserving its velocity components. 
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Schematically, this 1S done as follows: (The crystal regen- 
erates when the bullet passes plane E£E) 
A) Crystal Prior to Regeneration 


Direction of Travel 


Bullet The follwing planes 
ae are identical: 
i Ay- EB - YI 
B - F 
<= - G 
D-H 
A B Cc D 5B fe G HH I 


Disordered Crystal, 
Passed By Bullet 


B) Crystal During Regeneration 


a 
Plane E goes to 

F 1f f 

G ft tf 

H tf tf 

I tt tt 

Old planes A,B,C,D, 

A Ss © D 5 G H discarded 


Discarded 


HoaAw YP 


C) Crystal After Regeneration 
A p-——Bullet 


af 
a GC Pp -£ Be) 1 Planes F,G,H,I are newly 
Plane 
constructed as part of a 
perfect crystal lattice. 
Wa S. 
Plane I 


This scheme discards the portion of the crystal the bullet 
,has passed through and builds an undisturbed portion of the 


crystal in front of the bullet. 
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To build a crystal which will fit into the regeneration 
scheme the crystal generator must: 

(1) Construct.one plane at amtime sin sthe Y sdameciien sla 
places the first atom in the Y = O plane at the lowest value of 
X and Z (0,0) and continues placing atoms at increasing values 
of X for fixed Z, then increases the value of Z and repeats the 
process until the plane is completed. The next Y plane is con- 
structed in a similar manner and additional Y planes are 
constructed until the microcrystallite is completed. It assigns 
a number to each atom. A completed Y = O plane for the (100) 


Orientation has atoms positioned and numbered as follows 
ZZ @13 014 


e10 cal a 





Planes must be generated in this manner. 

(2) Build the following symmetry into the crystal. The left 
and right face are identical to each other and@identiical to one 
plane in between them. There is similar symmetry between the 
front and back faces as well as the bottom and top faces. 

This symmetry in the lattice allows the regenerator to 
Operate as follows: 

(1) It checks the location of each atom in the crystal to 


see 1f it is to be regenerated (shifted). 


a2 


(2) [if ae 1S toe be regenerated, a fixed number (IDX, IDY, 
or IDZ depending on the direction of regeneration) is subtracted 
from the atom number. This partially shifts the atom into its 
identical plane. 

(3) The regenerator then subtracts (or adds) DXT, DYT, or 
DZT from one of theatom's co-ordinates depending on the direction 
of regeneration. This completes the movement of the atom into 
its identical plane. 

(4) As previously mentioned the regenerator preserves the 
velocity components and potential energies of all atoms 
regenerated. 

(5) The bullet is regenerated ina similar (but not identical) 
Manner. 

In addition to being tailored to fit SUBROUTINE REGEN the 
lattice generator, through the use of a variable LCUT, breaks 
the lattice into an active core and a surrounding shell. This 
1s an absolute necessity due to the large size of the crystal 
which was dictated by the symmetry requirements of SUBROUTINE 
REGEN. (Orientation (100) contains 96 atoms, orientation (110) 
contains 60 atoms and orientation (111) contains 146 atoms). 
LCUT does this in the following way: 

1) Any atom in the active portion of the crystal is given 
an LCUT = O. This allows it to interact with any other atom 
gm the crystal. 

27 Any atom in shell is given an LCUT = 1 which allows it to 


interact only with core atoms and not with any shell atoms. 


B.S) 


The lattice generator subroutine also contains the constants 
XLL, XLF, ZLL, ZLF which are used to determine when the lattice 
1s to be regenerated. They are chosen to insure that the bullet 
always remains in the core and the lattice regenerates when the 
bullet 1s about to leave the core and enter the shell. 


2 
AS an example, a 20 Kev Xen y 


atom in the center of the (100) 
channel of silicon has a range of 491.7 LU and the program runs 
for 35 minutes and 39 secs. With the LCUT package inserted into 
the program the calculated range is 490.2 LU with a program 
running time of 13 minutes and 12 seconds. The LCUT package is 
Obviously effective for it conserves computer time without 
Sacrificing any accuracy in range calculations. 

In summary, the Lattice Generator 

(1) Contains all the constants necessary to regenerate the 
lattice, 

(2) It constructed, with appropriate tests, to insure that 
it fits into the scheme necessary for Subroutine REGEN, and 

(3) Contains the LCUT package which reduces the number of 
active atoms in the crystal, thus conserving running time. 

An alternative to building the lattice generator to fit 
subroutine REGEN would be to construct a smaller active crystal, 
eliminate the LCUT package and rewrite the REGEN subroutine to 
fit the lattice generator. This method was not chosen for the 
following reasons: 

(1) It was easier to modify the existing lattice generator 
to fit the old Subroutine REGEN than to start from scratch and 
build a new Regeneration Subroutine and Lattice Generator Sub- 


routine tailored for each other, and 
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(2) Simple modifications in the generator alone would not 
drastically alter the existing Channeling Program which had 
proved to be an accurate representation of the channeling pro- 


cess in fcc and bcc crystals. 
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APPENDIX B (CALCULATION OF FORCES AND POTENTIALS) 


A major change in procedure was introduced in the calculation 
forces and potentials. In previous models the interatomic po- 
tential function and force functions were of the exponential type, 

F = €xpm A + Bx) 
where A and B were empirically determined constants and x is 
the atomic separation. The atomic separation was determined, 
Substituted into the appropriate function and the force and 
potential determined. 

This computer simulation determined forces and potentials as 
follows. First, a set of potentials (in Rydbergs) derived from 


Interaction Potential Functions are read into the program for 


the lattice-lattice and bullet-lattice interactions. These 
Peventials Start at 0.0la,, (a, = Bohr radius = 0.5294) and are 
spaced every 0.Ola,, up to 5-0a,, interatomic separation. Secondly, 


a control card associated with each deck of potentials is read. 
This control card allows the extrapolation of a Born-Mayer type 
potential (if desired) from any point on the Interaction Potential 
Function Curve. Thirdly, the potential at the nearest neighbor 
distance is determined and then the potential is zeroed at that 
distance. All potentials at atomic separations smaller than the 
nearest neighbor distance are adjusted accordingly. All poten- 
tials at separations larger than the nearest neighbor distances 


were set equal to zero. This does not leave any potential free 
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regions in the lattice. Fourthly, potentials are converted to 
electron volts. Finally, the relation between the force and 
potential energy can be expressed simply as 
F(X) = - VU(X). 

Assuming that the potential is linear between any two successive 
points in the potential table, the force can be expressed simply 
as 

Bopees-—(Ulx) = UX + 4))- (const ) 
where the constant converts the force to proper units. 

The determination of a force or potential from the tables is 
just a simple matter of interpolation. For example, if the atomic 
separation between two atoms is 1.275a,, One would just make a 
linear interpolation between the values at 1.27a,, and 1.28a,,. 
These tables are used in both SUBROUTINE STEP which calculates 
forces and SUBROUTINE IONPE which calculates potential between 
the 10n and the atoms. 

Two points should be noted. First, the nearest neighbor 
separation is approximately 4.436a,,. Therefore, for a separation 
of 4.43h4a,, the interpolation is between the value in the table 
at 4.43a,, and the value at 4.4 36a. Secondly, the average force 
between two points is assumed to occur halfway between them. The 
interpolator for the forces accounts for this fact. 

One major problem remains unsolved in the construction of the 
force table (not in the interpolation of forces in SUBROUTINE STEP 
’ after the table has been constructed). While the potentials 
steadily increase as the interatomic distances decrease, the 


forces in the force table do not monotonically increase with 


o/ 


decreasing atomic separation. This is true only when potentials 
derived from the Interaction Potential Functions are used. No 
difficulty arises in the region where a Born-Mayer type of ex- 
ponential fit for the potentials is made. Here the forces behave 
as expected. 

The reason for this discrepancy is found in the calculation 
of the interaction potentials that are used by the computer pro- 
gram as data. The interaction potential is calculated exactly 
promo .0 lato 0.20a,, and every 0.20a,, to 5.0a,, (1.2 1On2 , One 
0.6a,, would be calculated exactly, 0.50a,, would not). Potentials 
in between these values are calculated by linear interpolation. 
Hence, the forces in the table would remain essentially constant 
in the region where the potentials were derived by linear inter- 
polation and would only change abruptly when a point is reached 
where the potential was calculated exactly. Thus while the 
overall force would increase as interatomic separation decreased 
there would be regions where the force stayed essentially con- 
stant or even decreased with decreasing atomic separation. 

Two main efforts were made to eliminate this discrepancy both 
of which began by considering only the potentials that were cal- 
culated exactly and discarding those that were determined by 
linear interpolation. First, subroutine ALI and ATSE from the 
scientific subroutine package of IBM were used to fit known po- 
tentials with a polynomial and using this polynomial, points in 
between the known potentials were interpolated. In the second 


method an exponential fit of the known potentials was made, and 
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an interpolator constructed as follows. Assume the potential P 


at two points xy and xX, are known and are to be joined by an 


exponential. We desire to know the potential at xX: 





KK =-X 
2 ; ; 

an = N where N 1s an integer 
el) PL = Aexp (-Bx, ) 
2) P. = Aexp(-B(x, +x, )) or - = Py exp (-Bx, ) 
Similarly 
(3) aa P, exp ( -BNAX) or In(P,/P, ) = -B*NAX) 
(4) Bax = - In(P,/P,))/N 


Hi 


Assume X, M-X Where Mis a constant O SM SN then 


Ee ate PL exp( -BAX:M) 


Both the polynomial and exponential fit of known potential 
yielded approximately the same results. In general, for atomic 


separations of less than 1.6a,, the forces increase with decreasing 


ay 


atomic separation as desired. For separations greater than 1.6a,, 
the forces act as expected, except that in the immediate vVacunaee 
(~.02a,,) of a calculated potential the forces decrease slightly 
and then increase again. Since there was little to choose between 
the two methods, the exponential interpolator was selected. The 
interpolator constructs the tables slightly faster (~ 2 secs) and 
requires less storage space. 

The effects of this problem on the force tables are as 
follows: 

a) te ach ERG effect. The forces behave as expected since 
the Born-Mayer approximation is constructed beginning at 1.2a 


pane nae 


b) Si =-Si (used as bullet-lattice interaction potential 


Hi 


in preliminary studies of (100) channeling). No effect, as a 
Born-Mayer approximation is constructed beginning at 1.2a). 
@)) Sul -si*! (used as lattice-lattice interactions for all 

directions of channeling). This force table contains the error 
discussed above for interatomic separations between 1.4 and 
3-0a,,- However, this error has negligible effect on calculated 
channeled ranges since (a) there is no error in the force table 
at atomic separations of 44a; the approximate nearest neighbor 
distance silicon atoms and (b) the bullet passes the lattice atoms 
SO quickly that when two lattice atoms have an interatomic spacing 
of less than 3-0a,, the bullet is no longer under their influence. 
While this error in the force table has virtually no effect 


on channeling range estimates, before the force table concept is 


used for other computer studies this problem must be solved. 
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APPENDIX C (PROGRAM PARAMETERS ) 


The following parameters have been added to the basic computer 
channeling model (Chan 69MODO). 
BPNT - Bullet-Lattice Potential Energy 
TPNT - Lattice-Lattice Potential Energy 
BFORCE - Bullet-Lattice Force 
TFORCE - Lattice-Lattice Force 
KA - an integer constant which identifies the point in the 
potential table where a Born-Mayer potential of form, 
he, begins 
BPCNST - the value of B in the Born-Mayer Potential approximation 
for the Bullet-Lattice Potentials 
TPCNST - the value of B in the Born-Mayer Potential approximation 
for the Lattice-Lattice Potential 
ICHECK, ICOUNT, NCOUNT - constants used in the exponential inter - 
polator between calculated interactions potentials in 
constructing the Potential Tables 
BDELX - a variable (BAX) (see appendix II) used in the exponential 
interpolator between interaction potentials in con- 
structing the Potential Table 
ITITLE - 80 alpha numeric characters used on a heading card for 
potential data deck 
LCUT - integer constant with value of either O or 1. Used to 
| determine which atoms are in the core (0) or shell (1) 


of the crystal 


Ay 


CVEDA - conversion factor used to convert electron volts per 
Bohr to newtons per meter 
CLU - conversion factor used to convert distances in LU to 


distances in units of 100a,, 


AFC - distance in 100a,, 


IFC - integer value of AFC 
CLU, AFC, and IFC are used in the following way. Assume 


interatomic spacing between atoms = 0.5LU 


Ane ==" (CrU\i oo DIST) 62 aon 
[PC = 256 
AMMAFC «of 4256. lMcorresponds to 2.5614a 
An IFC of 256 corresponds to 2.56a,, 

AFCI1 - nearest neighbor distance in units of 100a_, 

AFC2 - the reciprocal of the difference between the nearest 
neighbor distance and 443.0. Used in Subroutines STEP 
and IONPE. 

DFF1 - the difference in 100a,, between the atomic separations of 
the bullet and lattice atom and the nearest neighbor 
distance of silicon 

GETIME, SETIME - machine subroutines used to determine actual 
program running time. They terminate program execution and 
have cards punched before exceeding the running time allo- 


cated to the program by the computer. This allows a re- 


start capability. 
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Computer Program - SUBROUTINE DL100 


fe a ak kk KK ko ok SURROUTINE DOLLOO Ke KR ok tek a ke ak ok te de cK ok oo 
* THIS IS A LATTICE GENERATOR FOR THE OYJ AMOND LATTICE IN * 
x THE (100) CRICNTATICN, % 
* THE CPYSTAL NDEVELNPS IN THE ORDER X,FOLLOWED RY 7, * 
- FRELO@MPD FY Y¥ & 
Kok tokio tok kok kok kotokok koto tokok toto kok teak i kot koto ok kkk tet tok teak ake ak fc ak ok ke ak a ok 
En Semen be? Lacie ke Soh eCopt eran oe eee g LN 
COMMCN/COCM2/ROF ,ROF2 ~,ROFM,AC,PAC,PPTC,FTC,OFPYC, Porc, 
2 FM,PFIV, TOOT 
POON an tLe AER Aye decease 9 Chal EP 
COMMON/CC¥7/Q1 LSS, SPX ,S97,COX CHOY ,CN2 
SCOR ye ene! Ve xXe IOV eI NF 6200) -DXT YT, C7U,TECY wt BOY, 
FOMMON/COM22 /XLL »XLF»7LL»7LF 
ESS SELES OLE SESS SS CSCS CSS L OSES ELS ESS SS SS SS OLS SSE LE SELL LS Shee, 
~ THE FORESWIAG CONSTANTS GRE USER EITHFES BY THIS StR= * 
* RNUTINE sTO GENERATE THE LATTICE OR IN’ QTHEP SURREIITINES *® 
* fb REOUTPET . % 
sok KK tk kok kak tok gk kook gk kok akokoke ke tek teat ak ate tok ok ak ook ate ate ai ak kok ak de ak ak oe at ake ae oe ate ote 
XLL=C.¢S 
XLF=H=3.Cl 
Z7LL=C.S¢ 
7iF=3.01 
M=2 
L4=2 
TR=1 
1) X=] 
INY=4] 
TX=5 
TY=5 
| ees 
SC X=1.C 
SC vaoL eG 
S$C7=1.0 
Y=? ., C*xCC x 
NYV=2 6 C*xSC Y 
NZ T=2.0*SC7 
DY=0. > 
OX=0.5 
7 = 25 
KR KK tok i ik iotok gio totok tetototi tok iki tak kok ok oko skoak tet ote 2h ae tak oe ak ak ke fete aie ate ak ate 
* THE FOLLOWING SERTIFS CF NESTED AC LONOES GENERATE THE % 
* LATTICE AND INSURE IT CONTAINS THE ROQUIRED SYMMETRY * 
OR oq gk oko kok koko tk kok tok otoktok ok tok koko kok kook fa kokot kok kk tok ote teak ge ak ak ak 
JT=0 
Y=- SC Y 
0c AC J=L,TY 
Y=Y+SC Y 
KT=0 
7=-SCZ 
DN 50 K=1, 1.7 
7=7+SC7 
IT=0 
X=— SC xX 
DN 4CQ T=1,TX 
X= X+SC X 
LTTHIL T+ UT+KT 
IFC I TITHRACITT/2%*2)) 40530,40 
30 RXIM)=xX 
OYV(M)=Y 
R7(M)=7 
[D7(M)=5 
L=M 
M=M+}] 
40 |[T=[ T+] 
SCC KT=KT+4l 
Pete ltQO.5) GG TD 80 
HY 20C TA=LA4A oL 


43 


ee HHH 
% UI Hy 
+> + 
+r 
+r + 
HC + 
HT + 
+ + 
+e + 
Hc + 
+ + 
#+C + 
+p % on MK one 
+c + iT) q 
+ ca a a 
+ + tax a 
Hwy +t - FF 
+e ot —_ os 
Hy; Ht ey U 
Cie + > tt a (e 
CeO %* Ou + 
IAS HOT # ~ om 
- ¢ He ak Cc = 
Oo #u +* a ° 
Omg eC € al 
AIAN KHe- SH ° r 
e-¢ + mh + ee. 
Ce) H+ UC + wed om (" 
CoO Hus + ell e« 
CIN 5 ae ED - eae 
+ we CH a ol -—_ 
aN +H + <b BY 
OW Fn aw 3 Cn Or ae. 
SS + Wt eC) e — 
+t <> KR #a + [ae Sle 4 v) 
aml OO + Cw eC fo & 
a< ++ + %#O>* om © <I 
—uw <{ «<l <f + eee -~ 
>o  e HCO fre nc) = 
Yaoew— HOt eee AJ 
a SP all a | + ett rreu — 
xe Yo Il HTM e I Cre oe 
Ge WoW tie Hol + eC ee e CO 
UUimmmy, He SHOW YH mm mm, Ot 
e#H OC Ce wT at Ht Bayt i mer ee em 
CDC) pe he ee ee +t C Rew ww MZ 
ee¢ + + JU me _ ROH ~~ KP ALK ew Oo 
VI ede er et FHM S HZ HNRFaaaree LD 
wee wee eK HCI NN Ne HH De eZ OT 
wu X<PFnKRCeC Nace ye KFKOCCUUU C CRmWwz 
mmo oO ees WJ I OU HE OC OD es Lo os 
%* Wa 
Ne © Oct Tc Ft N 
Crea © ML Otkr OF 
OUAS OJ Sg + 
+t + 
% + tH 


44 


Computer Program - SUBROUTINE DL110 


KEKKER RRR RRR EH SURROUTENE SLL LO te Rt eR dR eK ee 
* THIS 'S @ tATTICE GENFRATOR FOR THE DIAMOND LATTICE IN * 
* THE (L190), CRISCNTATLONG ” 
> THe CRY 2*eh NEVebAP'sS IN THE CRDOER A, GELUGWER RY 725 * 


x Fepiirae 8 yY" ¥ * 
SSS SSE SLE LL SE SESE LES SELES SSL SS SSE SLES SES LL SS ES SSS SSS ES 2 2 2 5 33S f : 
COMMAN/OCMIT/RKX( 200) -2PY(200) «270 2090),LCUT (290),LL,LO 
COMMON IC CH2FRCE ,BMEZ » ROEM, AC , PL AC PETE, PTC DEPT ES EBTC, 
2 oar ve AS EU 71) T 
COMMAN /C ies PE ER ABS Sa CF, IDFEP 
COMMON/COMISRI LSS 4 SPX SP2 COX CHY (ama 
A ead IDX TAY, 107(290) gDXT,AYT, N7T,TOCK,T PNY, 
CIMMON/COM2P/XKLL »XLFE,Z7LL»7LF& 
‘See Le SS SP SEL SSE SES SPSS ES SSL SSS SSS SES 2S 2 2 2 SS SS 2S 2 2 2 St SSS SS SS 
* THE FOLLIWING CCNSTSNTS ARE WYSEN HF TTHER BY THIS SIR-= * 
* POUTINE TO GENFRATE THF LATTICE QR IN CTHFR SIJBROUT INES * 
Ee FECL ED, * 
SPSL SE SSS SSS ESSE EL EES SES SS SS 2S 2S SSS SSS SS SSS 2 SSS SS SSS SS SSS SSS 
XLL=C.4&C7 
XLF=2. 31 
rel =C 
i= sa 
M27 
(ox =" 
[ey=23 
I X=§& 
[ ¥=5 
I7=$ 
ONE 1.0/7S9R8T( 2.0) 
SC X=89 
S@y="f 
<@7=1.0 
NXT=2.,C*SC X 
OyYt=2.C*SGy 
N7T=2?,.C*SC7 
YX=P CR 
\F=0.0 
=a. * 
RR he aK ae oe 2c ae ae ae oe a ae ke ee oe ae ke he ae he oe oe ate he ee eae hea i ae ake ahead ke aie he atk ae ae ie ae ate ae ke ade ok ae XK oe oe 
* THE FOLLOWING SERTES GE NESTEN NO LONPS GENERATE THE x 
* LATTICE AND INSURE ITT CANTAINS THE REQUIRED SYMMETRY x 
SK WE whe he ae ae i oe ee Eafe ae a a ae ate a eof oe ee ake aoa eae ade ee ae ae ok tee ie ieee ote ie ic ake oe de ie ate ae i hea ok ake kk ok kk ote 
jT=C 
Y=- SC Y 
RP) eeo J=),1Y 
Y=Y+SC Y 
KR T=C 
7=~—SC 7 
oa Se K=)],[7 
Y=7+SC7 
I T=0 
X==- SC xX 
MA=M 
Y9 4C T=1,1X 
X=X+ SC xX 
IFC TT]$CIT/2)4*2) 21511,21 
LL TEC IT$€5IT/2) #2) 40,12 ,49 
12 TEC(KT—(KT/2)%*2) 40,30,40 
21 TFC IT-$€IT/2) 42) 22 ,40,22 
22 TFCKT=-(KT/2) #2) 30,40,30 
30 RX(“)=X 
RY(M)=Y 
R7(M)=7 
L =!4 
M=M+4] 
4Q [T=] T+}l 
MR=M—] 
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Computer Program - SUBROUTINE DLII11 


Kotoko Ro tk kkk = C SUIRROUTINE OLL1L1 See SSCS S Se Ste Se 2 2 eS 
* THIS 'S A LATTICE GENERATOR FOR THE OIAMOND LATTICE IN * 
* TRE CTL) CRIENTATICN. * 
* THE CPYSTAL DEVELOPS IN THE CRDER X,FOLLOWED BY 7, * 
* FOLLOWEN RY & 
SS SSSCE OCS CCS SCE CSE CS ES SSS SSS SSeS SS CSC TSS SSCS SS SS SSCS CLS Ss ee 
COMMON /CCOMIL/2X( 200) ~RY(200) AAP DOPE oT EC ORRT CURSE 
COMMCN/C CM2/&0F » ROE2 »ROEM, AC, PAC Cy,FTC,PFPTC,FPTC, 
FM,CFYV, TOOT 
COMMON /SCOMGSTX, LY ,gI7 oT XP TY Oe TZ PeSCX pSCY »SC7, IDEEP 
COMMON /COMT/R1 LSS 4 SPX 9 SOF yg CIXy COY CM 
eal IDX eINY,197(200) ,OXT,OYTsNZ7T -~FTPCX, TPO; 
COMMCN /COM22/XLL »XLF e7LL e7LF,YLLeYLF 
Se SCS SSCCSCCe* CT CSS ECS CL SSS SSC SS SST SSS SS SSS SS SS CSS SS TSS TSS Se TS 
* THE FOLLOWIAG CCNSTANTS ARF USED FITHER BY THIS SUB- * 
* POUTINE TC GENERATE THE LATTICE OR IN ATHER SIJBRCHT INES , 
* AS REQUTREN, 
ei RRR af Rak of att ak ak age ob aie ate a ate oe SF ERR ekl Ten te a vd. Fn MNT: ee Hee 
MLL=0.64°7 
ma F=Zz.131 
WEL=2.°5 
VYLF=4,4 
ZTLL=1.215 
ZLF=3,€F4 
M=? 
LA=2 
i= 1 
IN X=1 
INY=44 
M9 2 [=2,13C 
2 lua y j=5 
[X=5 
TY=7 
i = 
SC X=1-C/SQRT(2.0) 
SC Y=2.C/S9RT(3.0) 
Se7=SORTt lL. 9) 
Sc 7=°C 773.0 
NXT=2.C*SCX 
O¥T=2.C*#508T (330) 
O7T=2,0%*SC7 
DxX=0.0 
ap BRE 
=0.0 
PSPS SSL SSS SS SCSLSS LSS SSCS CC SLL SSS SSCL SSS SSeS OSCOSCS SS LEST es tS 
* THE FOLLOWING SERTES CF NESTED NO LNOFS GENERATE THE % 
* LATTICE ANP INSURE TT CONTAINS THE REQUIRED SYMMETRY x 
SES PLELE SS SSS CLS SSS SSS CS LSS SE LSS SSCL SC Se LC ST STS CLS SS St tt 2. 2s 
JT=0 
Y=- SC Y 
No 49 J=1L,IY 
Y=Y+SC Y 
KT=0 
JTS=5T# 9773 
7=-SC7 
90 50 K=l1, 
IF(G.£O0.©) G@ TR 8d 
CO CONTINUE 
Z=7+5C7 
IT=0 
X== SC Xx 
X=X+SC X 
IN=ITT#JTS4#KT 
TFC IN-€IN/2)%2) 40,30,40 
30 RX(M)=xX 
RY(M)=Y 
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Computer Program - Exponential Interpolator 


ao tok OE XPANENTIAL INTER PCLATOR Re RK 
Pen eT TPCE=LAVIITCE POTENTIALS AND DETERMINE 1° A BIRN—* 
* MAYER @PPRCXIMATION ITS NECESSARY “4 
eC OR OI IOI ok ta koi Gioia iio tok rotor ik tok iii togaioiotc kc kokog ca ae 
REAN(S, SO25VC TTI TLE CT) sf=1 20) 
READ(5,CA8IC)(TPNT(T) gl =1 920) 
REAN(5,SPLI5S)(TPNTCT) »f =21 4509) 
REAN(S5,°R82C) KA, TPCNST 
TF(KA) 405,405,410 
405 ICHF K=4éFC 
Cot 6 1¢ 
2 ROO OI RII IRI RRR IR RR RRR RI ROR ROR Ck to Xe Rete ate ak ae 


* CANSTRUCT A RORN=MAYER OCTENTIAL * 
-oPe ete SSSSLSLS SSL SESS SLES SSS SSS SSS SS SSE SSS 2 5 2 23 2:3 5 5325 2:5 5 595 | 
410 DO 41° 1[=KA,500 
415 TPNT(L)=TPNT(KA) #EXP( TPCNST#( (1/100 .)-(KL/10%~))) 
SoS SSS SLE SS SSS SSS SLL SLES S22 2:5 53-5 3 2 2:2 2-2-9 23 2 3 3 2:95 3 5 £2 3 2:3 5 2 SS Ss 
* MAKE AN EXECNENTISL INTERPOLATION AFTWEEN CALCULATED * 
* POTENTIALS 
S22 22 SSP SSS SSSSESS ES SSE SPS LES £2 S25 25 2:53 3 5 223 2 35 5 33:5 53 35:5 5 2:2 2 


ICHECK=KA 
416€ ICNOUNT=2C 
417 99 41° T=1,2C 
NC OUNT=TC CUNT#20 
BNELX=-ALCG(TPNT(NCQUNT) STENT (CIT COUNT ))/2C. 
41S TEONT(TCCUNT4+TI)=TONTO(TCOUNT) EXO ( (-1) *BDELX ) 
IFC ICHEC K=-NCOIUNT) 422 5,422,421 
421 ICNUNT=NCCUNT 
GN TN 417 
STC TCC TT SSCSCSOCSC§OCT COSC TST CCST CCT CST SS TTS SS TTS STL eS Te see eS 2 
* NETERMINE PCTENTIAL AT NEAREST NETGHBOR DISTANCE,7ZERNA * 
* POTENTIAL TABLE,AND CENVERT PATENTITAL TO ELFCTRAN VOLTS * 
CCCP CTCL SLCC LISS SLOT ISSCC TS ee TTS CLOTS SSS SCT TTS SLT TS es te te 
422 TIPNT=TENT(TEC) + (( TPNT(T ECHL) —T ONT (TFC) ®(AFC-TFC)) 
DO 420" B=1 . a0 
TPNT(T =O TONT(T )—TTPNT) ¥13.46 
IF( TONT(T ) 425,420,420 
425 TENT(T)=C.% 
420 CONTINUE 
ROAR OR ROR ROK kok kak kok 


222 222 555 2 22.29 2:2 2.9.2 23 25 2 23 23 2 2 2S 25 5 5 2 2 232 3 SS. 
x 


* AGONSTRUCER A LATTE E-LATTICE™FORCE YARLE 
ne aia aan" dlalaliie lll! is all siti!” “ainda 


iiss 
430 TEARCE( 4) ef TENT( 1) STONTII+1)P*CVEDA 
We Re whe 2 ode te ae ae oe ok he ode ote he ok oe ae ae ae ete ke oe oe ake oe ode ake oe ode ote otk dk te ate oe ate ok ok ate te otk ok ke ode 2 ie ate ok ok ok oe ok ok aie ote ke ok ok 
* READ RULLET-LATTICE PCTENTIALS AND NETERMINE TF A RORN- * 
* MAYER APDOROXTMATION TS NECESSARY a 
2 2 S222 SSS SS SSS SS SSS SSS SS SSE SSDS SSDS 2522.95 2 22.5.2. 2-2. 2:92: 2-2 S:9 2 2S: S| 
READ(5,CB25) (I TI TLEC(I) I =21 4,40) 
REAN(S,S5210) (BENT(T) g1=1 929) 
PEAD(5,CB1S)(BPNT(T) ,1=21 5509) 
RFAN(5,CR2C) KA BOECNST 
TS(Kh) 405,445,450 
445 ICHECK=46C 
5 eer ee 
Ce ae ean as ee 


* CONSTRUCT A BORN=-MAYER POTENTIAL 
ae, CT PND Se akg a i or ee er 
450 90 465 1=KA,500 
GEE BENT( 1) =BPNT(KA) #€ X2( BPCNST#((T/100~)-(KA/100.))) 
2 2 2 22 SF SSP S SPSS SSS SPSS SSS SPS SS SSS SES SSS SS SS SPSS SSS SS 2 S35 2 3.2 S | 
* MAKE AN EXPONENTIAL INTERPOLATION RETWEEN CALCULATEL o 
* POTENTIALS * 
(SPP SLL SSS SSS SLES S SESS ESS SSE SSS SSS SESE S SESS SES SSS S53 25 2 33-25 2°35 | 
ICHECK=KA 
45— ICOUNT=20 
457 N00 459 I=1,20 
NCOUNT=ICCUNT#+20 
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NELX=-A1 OSC BPNTENCOUNT) JAPNT(TCCUIT IIS 20. 
4590 RPNT(ICCUNT+I) =RONT(T COUNT) *EXC ( (=1) *BDELX ) 
TE( TCHFCK-=N\CCUNT) 462 74629461 
461 TCOUNT=AKC CUNT 
ie oh LET 


RICO TOK RR RR ROI ia Rokk tot kk tk ak katt ak te ak 
x OE TERM LE Ee Ten T ret BNE Res | NEIGHSOR OES TANCEsZERS % 
x POTENTTEL “SRL ,AND GENVERT PETENTHMt 1C ELEC” 2 GINEY Otis oe 
{OO RG OI II II IO IO RO I IO II I dai i ioi tok i ak kag kite a tek ak fc ak a 


462 RRPNT=RENTITEC) + ((BONT(TEC#L) -BPNT (TEC) )*(AFC-TEC) ) 
Nn 460 1=15 560 
RPNT(T)=(RONT(T)-BRPNT) #1326 
IF (RENT(T) 1475 9490 489 
475 AENT( Lae. °c 
LRO CONTINUE 
22 2 2S 22 SSS 22S PSS SSE SE SSS SSL ES SSE SS SE DSSS SESE LES ESS ES SS 2-2 22-25 
* CONSTRUCT A BULLET-LATTICE FORCE TABLE . 
8 ee erating, Theat acai ee eS Se ee ea ee 
GE5 AN 490 1=1,49 
Zee REGEGETTY =(90 ONTUT) -APAT(T#1)) #CVEDA 
eke ete ote kek ake ok ok ake ok ode ak ake a END OF INTERPNYATOO te ke ode te ak ake ake ake ake ae ke ake ate ak ok oe ake ote 3 
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Figure 4. The shape and sizes of the major 
channels in silicon. 
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Figure 5. 
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Note: 9 /em* = LOR = 15.5LU. 
(Ref 1,6F to. 5) 
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Figure 9. 
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RANGE PROFILE SECTION 
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Figure 14. 
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